%% ex4.m
% set dir: src/
close all;
clear; clc;

% Load image
imgC = imread("../data/cameraman.tif");
imgE = imread("../data/einstein.tif");

% convolution kernel
w = [-1 -1 -1 -1 -1;
    -1  2  2  2 -1;
    -1  2  8  2 -1;
    -1  2  2  2 -1;
    -1 -1 -1 -1 -1];

cReplicate = twodConv(imgC, w, 'replicate');
cZero = twodConv(imgC, w, 'zero');
eReplicate = twodConv(imgE, w, 'replicate');
eZero = twodConv(imgE, w, 'zero');

figure('Position',[100,100,900,300])
subplot(1,3,1)
imshow(imgC)
axis off
title('Origin')
subplot(1,3,2)
imshow(cReplicate)
axis off
title('Replicate')
subplot(1,3,3)
imshow(cZero)
axis off
title('Zero')
saveas(gcf,'../assets/cameraman_convolution.png')

figure('Position',[100,100,900,300])
subplot(1,3,1)
imshow(imgE)
axis off
title('Origin')
subplot(1,3,2)
imshow(eReplicate)
axis off
title('Replicate')
subplot(1,3,3)
imshow(eZero)
axis off
title('Zero')
saveas(gcf,'../assets/einstein_convolution.png')
%% twodConv
function g = twodConv(f,w,padding_option)
% Operate two dimensional convolution.
% - Input
% f(matrix): input grayscale matrix.
% w(matrix): convolution kernel.
% padding_option(string): padding method, providing 'replicate' and 'zero'
% option, default 'zero'

if nargin < 3
    padding_option = 'zero';
end

% determine size of kernel and image
[r, c] = size(f);
[rw, cw] = size(w);

% compute padding size
padH = floor((rw - 1) / 2);
padW = floor((cw - 1) / 2);

% padding boundary
fPadded = zeros(r + 2*padH, c + 2*padW);
switch padding_option
    case 'replicate'
        fPadded(padH + 1:end - padH, padW + 1:end - padW) = f;
        for i = 1:padH
            fPadded(i, padW + 1:end - padW) = f(1, :);
            fPadded(end - i + 1, padW + 1:end - padW) = f(end, :);
        end
        for i = 1:padW
            fPadded(padH + 1:end - padH, i) = f(:,1);
            fPadded(padH + 1:end - padH, end - i + 1) = f(:, end);
        end
        fPadded(1:padH, 1:padW) = f(1,1);
        fPadded(1:padH,end - padW + 1:end) = f(1, end);
        fPadded(end - padH + 1:end, 1:padW) = f(end, 1);
        fPadded(end - padH + 1:end, end - padW + 1:end) = f(end, end);

    case 'zero'
        fPadded(padH + 1:end - padH, padW + 1:end - padW) = f;
    otherwise
        error('Unsupported option, please check your input.')
end

g = zeros(r, c);
% conduct 2-D convolution.
for i = 1:r
    for j = 1:c
        region = fPadded(i:i + rw - 1, j:j + rw - 1);
        g(i,j) = sum(sum(region .* w));
    end
end
end
